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Abstract 



We study the dynamics of condensate formation in an inhomogeneous trapped 
Bose gas with a positive interatomic scattering length. We take into account 
both the nonequilibrium kinetics of the thermal cloud and the Hartree-Fock 
mean- field effects in the condensed and the noncondensed parts of the gas. 
Our growth equations are solved numerically by assuming that the thermal 
component behaves ergodically and that the condensate, treated within the 
Thomas- Fermi approximation, grows adiabatically. Our simulations are in 
good qualitative agreement with experiment, however important discrepancies 
concerning details of the growth behaviour remain. 
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I. INTRODUCTION 



The discovery of Bose-Einstein condensation (BEC) in trapped atomic vapors of 87 Rb |TJ, 
7 Li [0, and 23 Na || has initiated a period of intense experimental and theoretical activity. 
A great deal of information is now available about the equilibrium properties of these novel 
systems [|J, but much remains to be understood about their nonequilibrium behavior. One 
of the most basic aspects concerns the nonequilibrium growth of the condensate which 
occurs in the process of cooling a nondegenerate trapped Bose gas to a final temperature 
below the BEC transition. This important problem was addressed even before the first 
observation of BEC in trapped atomic gases and has interesting implications for the 

general problem of second-order phase transitions, from superfluidity in liquid 4 He || to 
problems in cosmology M. 

Up to now, the most detailed study of condensate formation was carried out using a 
gas of 23 Na atoms confined within a highly anisotropic cigar-shaped trap |1U||. In these 



experiments, the sodium atoms were evaporatively cooled to a temperature just above the 
critical temperature and subsequently quenched by applying a rapid rf sweep. The latter step 
removes all, or at least a large fraction, of the atoms above a certain energy, after which the 
Bose gas relaxes to a new equilibrium state below the critical temperature. The growth of the 
condensate during the equilibration process was monitored using a nondestructive imaging 
technique which provided a direct measure of the size of the condensate as a function of 
time. In this way, the characteristic time scale for the growth of the condensate could be 
determined, and for the particular system studied, was found to be of the order of 100 ms. 

A theoretical description of these experiments requires a theory that can account for the 
coupled nonequilibrium dynamics of both the noncondensed and condensed components of 
a trapped Bose gas, and includes in particular the collisional processes which transfer atoms 
between the two components. Thus far several such theories have been developed, which 
roughly speaking fall into two categories. One class of theories focuses on describing the 
dynamics of the average value of the order parameter for BEC, i.e., the condensate wave 
function, whereas the other incorporates also the fluctuations around this mean value. The 
latter of course, becomes important when the fluctuations are large compared to the mean 
value, i.e., close to the critical temperature. This is analogous to the situation in laser 
theory JTTJ. 

A theory that describes both the average value for the order parameter as well its fluc- 
tuations can be obtained in two, essentially equivalent ways. First, one can start from a 
master equation for the many-body density matrix and derive an equation of motion for the 
one-particle density matrix by means of a perturbative treatment of the interactions. This 
was the route followed by Gardiner and Zoller [O, in a series of papers. Second, one can use 
field-theoretic methods to obtain a nonperturbative Fokker-Planck equation that describes 
the nonequilibrium dynamics of the gas. This was the formulation developed by Stoof [JT~3 



These two approaches in principle yield a description of the nonequilibrium dynamics that 
is capable of obtaining the complete probability distribution for the order parameter. 

Alternatively, a theory describing the dynamics of the mean-field value for the macro- 
scopic wave function can also be obtained in a straightforward decoupling approach, which 
has been implemented by Kirkpatrick and Dorfmann [H|, Proukakis et al. ||15|| , Walser et 



al. [16], and in most detail for the trapped case by Zaremba, Nikuni, and Griffin [17|. In 
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this approach, one assumes the order parameter to be nonzero at all temperatures, and 
decouples the hierarchy of equations of motion that exists for the correlation functions of 
the second-quantized field operators. Thus one again obtains a perturbative expansion for 
these equations of motion. 

The first quantitative calculations of condensate growth for trapped Bose gases were car- 
ried out by Gardiner et al. |Tj|, and although good qualitative agreement with experiment 
was found, a number of quantitative discrepancies remained. For example, the reported ex- 
perimental growth rates were up to a factor 30 larger than the initial theoretical results, and 
had a temperature dependence opposite to that predicted flT0|. By removing some simpli- 



fying approximations in subsequent calculations, the theoretical results were improved, but 
discrepancies of up to a factor of 3 still remained in some cases. From a purely theoretical 
point of view, one can attribute some of these discrepancies to the approximations made in 
the calculations. First, the dynamics of the noncondensate was to a large extent neglected. 
Although the time evolution of the occupancy of low-lying states was included in the simu- 
lations, the high energy states were represented by an equilibrium particle reservoir having 
a fixed chemical potential. This latter assumption is inconsistent with the nonequilibrium 
initial state established by the experimental quench procedure. Second, the effect of the 
mean field of the condensate on the noncondensate was included rather crudely by a linear 
rescaling of the low-energy density of states of the noncondensed atoms. 

Our aim in the present paper is to improve on these calculations by taking fully into 
account the relaxational dynamics of the thermal, or noncondensed, component that takes 
place in the presence of the mean field of the condensate. We do this by starting from the 
above mentioned theories describing the growth process, which provide us with a nonlinear 
Schrodinger equation for the condensate and a kinetic equation for the thermal compo- 
nent. This coupled set of equations is still difficult to deal with and a number of phys- 
ically motivated approximations are made to simplify the problem. We assume that the 
condensate grows adiabatically, having an equilibrium spatial distribution determined by 
the instantaneous number of atoms in the condensate. This assumption is also made in 
earlier work W^. The noncondensate is treated by solving a semiclassical Boltzmann equa- 



tion |13|,|T^,|T^,|T9| in the ergodic approximation, which again has been used previously by 
numerous authors |||20|-^]. These assumptions allow us to obtain numerically a detailed 
description of the growth of the condensate, including the effects of both the dynamics of 
the thermal cloud and its mean-field interaction with the condensate. 

The paper is organized as follows. In Sec. |TI| we summarize the theory of the nonequi- 
librium dynamics of a trapped Bose gas as developed previously. In Sec. |II| we introduce 
the central assumption, the ergodic approximation, that allows us to numerically solve the 
Boltzmann equation. In addition, we briefly discuss the adiabatic approximation for the 
condensate. In Sec. |V| we treat in some detail particle number and energy conservation. 
Sec. [V] introduces the Thomas-Fermi approximation and gives some analytical results for 
the density of states and other quantities of interest. The numerical solution of our kinetic 
equations is discussed in Sec. |VT| and our results for the growth of a condensate are presented. 
We end in Sec. IVIII with a discussion and an outlook. 
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II. NONEQUILIBRIUM DYNAMICS 



As discussed in the previous section, the nonequilibrium dynamics of a trapped Bose 
gas is governed by a set of equations for the condensate and noncondensate components. 



These equations have been presented in various forms |l^^3| -[T7|1 , but they all describe 
the coherent dynamics of the gas due to mean-field interactions, as well as the incoherent 
dynamics associated with atomic collisions. The equations and notation we use are taken 
from Refs. fll3| , |TT . 



The noncondensate is treated using a semiclassical Boltzmann equation for the phase 
space distribution function /(r, p, t). This semiclassical description is justified when the 
largest level spacing in the external trapping potential is small compared to the thermal exci- 
tation energy. Moreover, mean-field interactions are included at the level of the Hartree-Fock 
approximation. In this situation, the quantum kinetic equation for the thermal excitations 
takes the form (EJ|I7| 

df(r d ^ t] + £ ■ V/(r, P, t) - W(r, t) ■ V p /(r, p, t) = C 12 [f) + C 22 [f) . (1) 

Here, the effective potential U (r, t) = U ext (r) +2g[n c {r, t)+n(r, t)] is the sum of the external 
trapping potential U ext and the self-consistent Hartree-Fock mean field. The latter is de- 
termined by the condensate density n c (r, t), defined below, and the noncondensate density 
h(r, t) given by 

fl ( r - t ' = /(W /(r ' p ' i) ' (2) 

As usual, we treat the interactions in the s-wave approximation which results in the bare 
interaction being replaced by a contact interaction with an effective coupling constant g = 
Aixh 2 a/m proportional to the s-wave scattering length a. The effective coupling constant is 
in fact equal to the two-body T-matrix, and to emphasize this connection, it is denoted by 
T 2B in some works fl3|]. The collision terms appearing in Eq. ([!]) are given by 

4vr 2 r dp 2 f dp 3 r dp 4 3 

x6(E + E 2 -E 3 - E 4 ) [(1 + /)(! + f 2 )hh ~ // 2 (1 + / 3 )(1 + h)\ , (3) 



n \f\ An 2 [ dp 2 f dp 3 f dp A 3 

CM - T 9n c J j^J j^J j^(2Kh) 5(mv c + P2 - P3 - P4 ) 
x6(E c + E 2 -E 3 - E 4 )(27ihf[5(p - p 2 ) - 5(p - p 3 ) - 5(p - p 4 )] 

x[(l + / 2 )/3/4-/ 2 (l + /3)(l + /4)] , (4) 

with / = /(r, p, £), and fi = /(r, p i; t). We note that Eq. @) takes into account the fact 
that a condensate atom locally has an energy 

E c (r, t) = /i c (r, t) + imv c 2 (r, t) , (5) 
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a momentum mv c (r,t), and a chemical potential /i c (r, t). These quantities are defined 
explicitly below. In addition, the energy of a noncondensate atom in the Hartee-Fock ap- 
proximation is 

E(r,p,t) = ^ + U(r,t). (6) 

The energy variables Ei appearing in Eqs. @ and (H) are defined as Ei = E(r,Pi,t). 

In contrast to the thermal cloud, the dynamics of the condensate is determined by a 



time-dependent dissipative nonlinear Schrodinger equation [13,17 



ih ^A = j-^ + f/ cxt (r) + g [2n(v,t) + n c (r,t)} - iR(r,t)\ *(r,t) , (7) 
where the dissipative term, i.e., R(r,t), is given by 



g 2 



(2tt) 5 ^ 



/ Yl dpi5(mv c + p 2 - p 3 - p 4 ) 



i=l,4 



x5(E c + E 2 -E 3 - E,)[5( Pl - p 2 ) - 5( Pl - p 3 ) - 5( Pl - p 4 )] 

x[(l + / a )/3/4-/ a (l + /a)(l+/4)]. (8) 

The appearance of this dissipative term in Eq. (0) is a consequence of the collisional pro- 
cesses, described by Cu, which have the effect of transferring particles between the conden- 
sate and noncondensate. The dissipative term is needed in order to ensure overall particle 
number conservation of the entire system. At a more fundamental level, the condensate 
wave function is determined by taking the expectation value of a Bose field with respect to a 
probability distribution that satisfies the Fokker-Planck equation mentioned previously . 



The Langevin equation one derives within this formulation has the form of a dissipative 
nonlinear Schrodinger equation with a noise term. It only reduces to Eq. ([7|) in a mean-field 
approximation. This points to the need for exercising care in interpreting the order parame- 
ter occuring in the nonlinear Schrodinger equation as the condensate wave function. Due to 
the underlying £7(1) gauge invariance associated with strict particle number conservation, 
the expectation value is in fact always equal to zero as a result of the diffusion of the global 
phase of the condensate wave function. In first instance this effect can be neglected and we 
are then effectively treating the system as if the U(l) gauge invariance is explicitly broken. 
It is convenient to rewrite Eq. (|7|) in terms of amplitude and phase variables defined by 



\I/(r, t) = J n c (r, t) exp[i6(r, t)]. Substituting this form of the wave function into Eq. (0), we 
obtain 

' : '""' r ' n + V[v c (r,tKM)] = ~R(v,t)n c (r,t) , (9) 



and 



m — h V 

at 



mv c (r,t) 2 



. (10) 
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Here, we have defined the local chemical potential and superfluid velocity by 



h 2 V\ n c (r,t) 

fJL c {T,t) = U cxt (r)+g[n c (r,t) + 2h(r,t)] - )L , (11) 

2m Jn c (r,t) 



and 



v c (r,t) = ^V0(r,*), (12) 

respectively. It can easily be shown that this set of equations for the condensate and thermal 
cloud is consistent with the conservation of the total number of particles in the system W? 



III. ERGODIC APPROXIMATION 

Our main objective in this paper is to apply the kinetic theory formulated above to 
the problem of condensate formation. In order to make progress we introduce a number 
of additional approximations. The first and most essential, is the assumption of ergodic- 
ity |6|j20|-p2l which has been widely used in the literature on kinetic theory. This assumes 
that equilibration of atoms within one energy level occurs on a much shorter time scale than 
equilibration of atoms between different energy levels. With this assumption, all points 
in phase space having the same energy are equally probable, and the distribution function 
therefore only depends on the phase space variables through the energy variable E(r,p,t), 
i.e., /(r,p,t) = g(E(r, p, t), t). In equilibrium this is certainly correct, but the assumption 
requires justification for any particular nonequilibrium application. Unfortunately, we are 
not aware of any explicit checks that have been made which might indicate that the as- 
sumption is correct for the situations we wish to consider. Nevertheless, it appears to be 
physically reasonable that for quantities that vary on a time scale of the order of several 
collision times, the approximation is sufficiently accurate. 

The ergodic approximation allows us to derive a simplified kinetic equation for the energy 
distribution function g(e,t). This is accomplished by means of the relation 

p(e, t)g(e, t) = j - E(r, p, f))/(r, p, t) , (13) 

which shows that the phase-space projection defined on the right-hand side yields the product 
of g(e,t) and the density of states 

" (e - i) = /(W a(£ - £(r - p - t)) 

™3/2 n 

= -j= ? / drJe - U(t, t) . (14) 

We note that the density of states is defined on the variable energy range U min (t) < e < oo 
where U m i n (t) is the minimum value of U(r,t) at time t. The time dependence of the 
density of states is one of the aspects distinguishing the present development from previous 
work HIO]. 
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We now apply the phase- space projection to the kinetic equation in Eq. ([[]). As a result 
of this operation, the streaming terms in the Boltzmann equation, i.e., the second and 
third terms on the left-hand side of Eq. ([I]), cancel each other. Only the projection of the 
time-derivative term survives. This results in 

f drd P x( cv ,^ d f( r >P,t) ( ,^(e,t) A dg(e,t) 

J ww 6{€ - E(r ' p ' t)] ^r- = p(e ' t] ^r + Pw(e ' t] ^r ■ (15) 

Here, we have introduced a weighted density of states 

drdp dU(r,t) 



(2Tihf v v " dt 
m 3 / 2 f I — ~dU{r,t) 



V2Tr 2 h 3 Ju< 

This quantity depends explicitly on the time derivative of the noncondensate potential which 
in turn is determined by the time derivatives of both the condensate and noncondensate 
densities. Some formal details regarding its evaluation are given in the Appendix. 
Noting that 



dp(e,t) _ dp w (e,t) 
dt ~ de 

Eq. flTSP can be written as 

drdp <9/(r, p, t) d . . d 



(17) 



f drdp C/ r,p,t 

J (2# 6(6 - E(r ' P ' fiT" = Ft {P9) + We M ■ (18) 



(27rfl) 

We thus arrive at the projected kinetic equation 



dt ^ + We (Pw ^ = /l2 + 122 ' (19) 



where the phase space projections of the collision integrals are defined as 

/l2M) = / ll^^^P'^^t/] ( 20a ) 



'221 



t)=J 7^<K e - E ^ P' *))C*\f\ • (20b) 



The result in Eq. (|19D is the kinetic equation that we solve numerically. 

We now derive in some detail explicit expressions for the collision integrals in Eq. (pO"D- 
Although an expression for J 2 2 was given in earlier work |2(J , we present here an alternative 



derivation which can also be adapted to the case of the I\2 collision integral. For the I22 
collision integral we have 
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Iw(ei,t) = T^p^To / de 2 J de 3 J de^S^ + e 2 - e 3 - e 4 ) 
x K 1 + + 52)53^4 - gi92{l + #0(1 + 94)] 

X / rfr f II / rf P^ 5 (Pl + P2 - P3 - P4) 



x«y(ei - E 1 )5(e 2 - E 2 )5(e 3 - E 3 )5(e 4 - E A ) , (21) 

where we have introduced the short-hand notation gi = g(ei,t). We consider first the 
momentum integrals in Eq. ( |21| ) which, with the replacement P3 — > — P3 and p 4 — > — P4, 
can be written as 

J 22 = J dpi J dp 2 J dp 3 J dp 4 5(px + p 2 + P3 + P4) 

X<f(ei - E 1 )5{e 2 - E 2 )6(e 3 - E 3 )5(e 4 - E 4 ) 

=/(^n 4 /^ e ^-^- ^ 

In obtaining this expression, we have introduced a Fourier representation of the momen- 
tum conserving delta function. Performing the integrals in Eq. (^) with respect to the 
momentum variables, we obtain 

t (a \ifriQt tt\\ f d € s in £Pisin£P2sin£p 3 sin£p 4 
J 22 = (47rm) ^ ]Q 0(e< -U)j J - 4 , (23) 



\i=l 



where now it is understood that pi = y 2m(e, — U). The product of theta functions can be 
replace by 9(e m - m — U), with e min the minimum value of the four energy variables. Performing 
the remaining integral with respect to the £ variable, we find 

J 22 = (27r) 3 m 4 0(e min - U) |pi - p 2 + p 3 +p±\ - \pi - p 2 + P3 - Pa\ 

+ \pi-P2-P3-Pi\- \Pi -P2-P3 + Pi\ 

+ \P1+P2+P3-P4~ \Pl +P2-P3~P4\ 

+ \P1 +P2 ~P3 +Pi\ - \Pl +P2 +P3 +Pa\ ■ (24) 

This expression is valid for arbitrary values of the momenta but simplifies when energy 
conservation is taken into account. Since the energy conserving delta function 5(e\ + e 2 — 
£3 — q) in Eq- ( |21| ) imposes the constraint p\ + p\ = p\ + pi, Eq. (pi[) can be reduced to 



J 22 = 4(27r) 3 m 4 fl(e min - U) ^/2m{e min - U) . (25) 
Substituting this expression for J 2 2 into Eq. (|21|), we finally obtain 



,3„2 



h^iit) = ^3^7 J de 2 J de 3 J de 4 p(e min )5(ei + e 2 - e 3 - e 4 ) 

x [(1 + g x ){l + g 2 )g 3 g A - 9l g 2 (l + g 3 )(l + g 4 )) , (26) 



where we have used the definition of the density of states in Eq. (|14[) . This is precisely the 
result obtained by Snoke and Wolfe |§ and Luiten, Reynolds, and Walraven pOf , using a 
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different method. We note that if all energies are expressed in units of hu, the I22 integral has 



an overall factor of (a/l) 2 cu where I = yh/mtu is the average harmonic oscillator length. This 
factor defines a characteristic time which can be used as the time unit in the simulations. 

The J12 collision integral can be dealt with in a similar way if the superfluid velocity v c 
in Eq. (§) is set to zero. The validity of this approximation follows from our assumption 
that the condensate grows adiabatically. The magnitude of the superfluid velocity v c is then 
typically of the order of R(t), where R(t) is the radius of the condensate. This velocity 



is small compared to the characteristic velocities p/m « y2/csT '/m of the thermal atoms 
participating in a collision, which justifies the neglect of mv c in Eq. (^). The expression for 
/12 then reads 

112(^1, t) = 7^6^7 / drn c (r,t) J de 2 J de 3 J de 4 6(E c (r, t) + e 2 - e 3 - e 4 ) 

x [5(e 1 - e 2 ) - 5(e 1 - e 3 ) - 5{e x - e 4 )] [(1 + g % )g 3 g± ~ #2(1 + 9z)0- + #4)] 

x J dp 2 J dp 3 J dp 4 5(p 2 - p 3 - p 4 )5(e 2 - £ 2 )5(e 3 - E 3 )6(e 4 - E 4 ) . (27) 

If we now define J i2 analogously to J 22 , we have 

J12 = J dp 2 J dp 3 J dp 4 5(p 2 + P3 + P4)5(e 2 - £ 2 )5(e 3 - E 3 )5{e± - E±) 



d£ sin £p 2 sin £p 3 sin £p 4 



= (47rm)^(e min - 17) } ^ ^ 
= 87T 2 m 3 0(e min - U)S(p 2 ,p 3 ,P±) , (28) 
where e min is the minimum value of e 2 , 63, and e 4 , and 

S(P2, P3, Pi) = ^ [ s g n G°2 + P3 - P4) + Sgn(p 2 - P3 + P4) 

-sgn(p 2 + p 3 + p 4 ) - sgn(p 2 -P3- P4)] ■ (29) 

Note that this is a boolean function which takes on values of and 1. Inserting the expression 
for J\2 into Eq. ([27D, we finally obtain for the 7 42 collision integral the result 

7i 2 (ei,£) = 3 ^ 7 y de 2 y de 3 J de 4 [5{e 1 - e 2 ) - <J(ei - e 3 ) - (J(e a - e 4 )] 
x [(1 + ^2)^4 -^2(1 + £3) (1 + 04)] 

x/ drn c (r,t)S(p2,p 3 ,p 4 ) 5{E c (r,t) + e 2 - e 3 - e 4 ) . (30) 

JU<e m ;„ 



A comparison of this expression with 7 2 2 in Eq. fl26|) shows that the remaining spatial integral 
acts as an effective density of states for scattering into the condensate. It can be evaluated 
analytically in the Thomas-Fermi approximation for the condensate, as shown in Sec. [V]. 
The kinetic equation in Eq. (|19D and the projected collision integrals in Eqs. (|26|) and ( |30"D 
are the main results of this section. 
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Before closing this section we point out one difficulty encountered when a numerical 
solution of Eq. fll9D is attempted. As discussed following Eq. (|14 ), the time dependence of 



the mean field potential U(r,t) implies that the density of states in Eq. (|H|), and hence the 
energy distribution function g(e, t), are defined on a variable energy range. To eliminate this 
variation, it is convenient to introduce the shifted energy variable 

e = e - U min (t) , (31) 

which leads to a fixed energy range < e < oo. The density of states in terms of this new 
energy variable is given by p(e,t) = p(e + U m m , t) = p(e,t). With e and t as independent 
variables, the kinetic equation in Eq. (|i~9|) , can be rewritten as 

Here, both p(e,t) and p w (e, t) are defined by making the replacement 

U(r, t) -> U(r, t) - U min (t) = U(r, t) , (33) 

in Eqs. fl!H| ) and (0). Similarly, from the definition of the collision integrals in Eqs. (|^) and 
(|), it can also be seen that the change of energy variable leads to the replacement of U(r, t) 
by U(r,t) in this case as well. Thus, the final kinetic equation in terms of the e variable 
is unchanged in form from the original equation. We will henceforth drop the overbar on 
the functions defined in terms of e, with the understanding that the shifted potential U (r, t) 
is to be used wherever the potential appears in the original expressions. The possibility of 
using a fixed energy range in the solution of the kinetic equation simplifies the numerical 
calculations considerably. 



IV. COLLISIONAL INVARIANTS 



In this section we explicitly consider two important quantities that should be conserved 
as the Bose gas condenses and equilibrates, namely, the total number of particles and the 
total energy of the trapped Bose gas. Together, they determine the final equilibrium state 
of the Bose-condensed gas, i.e., the number of particles in the condensate, its chemical 
potential, and the temperature of the vapor. 



A. Particle Number Conservation 

The time rate of change of the total number of particles consists of the time rate of 
change of the number of condensed particles plus the time rate of change of the number of 
noncondensed particles. Because the number of noncondensed particles is given by N(t) = 
(2nh)~ 3 J drdpf(p,T,t) = J dep(e)g(e), the time rate of change of N(t) can be found by 
integrating Eq. (|T9| ) over energy. We thus find that 

^ = /*MM), (34) 
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where it is easily checked from Eq. (|26"D that 



Jdel 22 (e,t) = 0. (35) 

Note that we have assumed here that hm e ^u mitl Pw( e )<?( e ) = 0. A finite limiting value can 
arise if g(e) approaches an equilibrium Bose distribution with a chemical potential fi = U m i n 
at long times, together with a weighted density of states which depends linearly on e — {7 min 
for energies close to U min . However, at any finite time in the growth process, it is safe to use 
the zero limiting value. This is always the case when the equilibrium chemical potential lies 
below U m i n . 

To get the time rate of change of the total number of condensate particles, we integrate 
the continuity equation, Eq. (d), over space to find 

dN c 2 r 

= —frJ drR(r,t)n c (r,t) 

= - J del l2 (e,t). (36) 
Combining this with Eq. ( pl[ ) leads to 

which demonstrates that the total number of particles is indeed conserved. 



B. Energy Conservation 



We now consider the conservation of the total energy of the system. The total energy is 
given by 



E, 



tot 



drdp ( p 



(2tt^) 3 
+ J dv^*(r,t) 



^ + U cxt (r) + g [n(r, t) + 2n c (r, t)} \ f(r, p, t) 



h ^ T T / \ 9 , X 

~2— + C/ext(r) + 2^(r,t) 



*(r,t) 



(38) 



The first term is the semi-classical expression for the total energy of the noncondensate. It 
contains the kinetic and external potential energy, and the Hartree-Fock mean-field interac- 
tion energy of the noncondensed cloud interacting with itself and with the condensate. We 
note that the self-interaction term is reduced by a factor of two relative to the condensate 
term to avoid double counting this contribution. 

The second term in Eq. (|38|) is the total energy of the condensate which contains the 
wave function \l/(r, t) with normalization 



J dr\V(r,t)\ 2 = N c (t). 



(39) 



It consists of the kinetic energy, the potential energy, and the mean-field energy due to 
the interaction of the condensate with itself. The mean-field interaction of the condensate 
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with the noncondensate has already been included in the expression for the energy of the 
noncondensed cloud. We now show that this total energy is indeed conserved during the 
growth process. 

Taking the time derivative of Eq. fl38|) leads to the following expression, 
dEufi f drdp f p 2 \df(r,p,t) 

+ J dr ° ^' > + C/ ex t(r) + 5 k(r, t) + 2n(r, t)] | ¥(r, t) 

+ y dr **(r, t) + U ext (r) + 9 [n c (r, t) + 2n(v, t)] j — ^ (40) 

The first term in Eq. ( f4"0"|) can be rewritten as 

/" rfrrf P P / .^df(r,p,t) f drdp rfl_i_r» rm 

7 (2^)* s( r ' p ' * )_ aT- = y (2^ E(r ' p ' t} (^ + ^1/1) 

= h E ^I(^^ (41) 

where, to obtain this result, we have used the kinetic equation, Eq. ([!]), and the fact that 
the C22 collision integral conserves energy. 

If we again assume that the condensate grows adiabatically as atoms are fed into it from 
the noncondensate, the condensate wavefunction ^(r, t) is a solution of the instantaneous 
Gross- Pit aevskii equation 

f h 2 V 2 1 

+ U eXt + g[n c (r,t) + 2n(r,t)]}*(r,t) = E c (t)V(r,t), (42) 



2m 



with a time-dependent energy eigenvalue E c (t). For this spatially independent condensate 
energy, Eq. (|4]) reduces to 

drdp <9/(r,p,t) gjfo) 

#(r,p,t) 7T7 = £ c (t)^— . (43) 



(2tt7z) 3 v y 0i cv 7 ^ 

Inserting this result and Eq. (p|) into Eq. (|40[) , the latter is easily seen to yield 

= , (44) 

due to the conservation of total particle number. Thus the assumption of adiabaticity is 
sufficient to ensure that the total energy is conserved. However, one can also show the con- 
servation of energy exactly, without assuming adiabaticity, by making use of the dissipative 
nonlinear Schrodinger equation in Eq. (|7|). 
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V. THOMAS-FERMI APPROXIMATION 



The assumption that the condensate grows adiabatically implies that the dynamics of 
the condensate itself is being neglected, apart from its trivial time dependent normalization. 
In particular, we are ignoring the possible excitation of internal collective oscillations. How- 
ever, at the temperatures of interest in the growth process, these excitations are strongly 
damped and we expect the condensate to remain in a relatively quiescent state which is 
well approximated by the quasi-equilibrium solution of the GP equation. Indeed, in the 
experiments there is no evidence of condensate oscillations, although the thermal cloud has 
been observed to oscillate at twice the harmonic oscillator frequency of the trap in some 
cases. 

For a large number of condensate atoms, a good approximation to the equilibrium wave 
function is provided by the Thomas-Fermi approximation which neglects the kinetic energy 
in the GP equation. In this situation, the condensate density is given by 

n c (r, t) = - [fi c {t) - C/ ext (r) - 2gn{r, t)] . (45) 

Of course, this expression is only valid if the right-hand side is larger than zero; otherwise, 
n c (r,t) = 0. The last term on the right-hand side reflects the mean-field interaction of 
the condensate with the thermal cloud. Since the latter has a small density relative to the 
condensate, its effect on the spatial distribution of the condensate is small (2gh <C /i c ) and 
we therefore neglect it when determining the condensate density. By the same token, we 
shall neglect the mean-field interaction of the noncondensate with itself. Strictly speaking, 
these approximations lead to a violation of total energy conservation, but the error will be 
very small since the bulk of the mean-field energy, which resides within the condensate itself, 
is still taken into account. In principle, these contributions can be included in our treatment 
as shown explicitly in the Appendix. However, because these corrections are small, we have 
decided to neglect them in our numerical calculations. 

It should be noted that the Thomas-Fermi approximation is to some extent dictated 
by our semi-classical treatment of the noncondensate atoms, since it avoids a potential 
problem associated with the placement of the condensate chemical potential fj, c relative to 
the minimum energy available to the thermal atoms, i.e., U m i n = min[?7 ex t + 2g(n + n c )]. For 
small condensate densities, it is possible that the GP eigenvalue \x c lies above this minimum 
value which is clearly impossible if a full quantum treatment of the excited states is retained. 
In the Thomas-Fermi approximation there is no such problem since the chemical potential 
is exactly equal to U min . 

Given these approximations, the time-dependent condensate density profile becomes 

ri c (r,*) = -M0-£/ext(r)] , (46) 
9 

where the external potential is taken to be a general anisotropic harmonic confining potential, 
U e xt( r ) — J2i mu! i r i /2. This expression for the density is again only meaningful when 
^ext( r ) < Hc{t)- The chemical potential of the condensate is given by 



/W) = — 



15W (*)y 



2/5 

(47) 
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where ui = (ui U2UJ3) 1//3 and I = yh/mu. The potential experienced by the thermal atoms 
is then 



U(r,t) 



(2fi c (t)-U cxt (r), ifn c ^0 
\U ext (r), ifn c = 0. 



(48) 



The minimum value of this potential is // c (t) and occurs on the boundary of the condensate. 

Three additional important quantities can also be calculated analytically. The first two 
are the density of states, and the weighted density of states, i.e., p(e) and p w (e) respectively. 
For the former we get 



in 



3/2 



drJe-U(T,t) 



■kTilo 



dyy 2 e - 



I2(e-fi c 



+ y 2 



+ 



dyy 2 6 y 



2/i c \ / 2(e + // c ) 2 

ho J V ftfi 



TlTlUO 



-[/_(e)+/ + (e)] 



(49) 



The integrals /-(e) and i+(e) are standard, and are given by 

— log (a: + u 



v?_x a^u_x a 2 



f- e 



Ue) 



4 



u%x a + u + x a\ 

H 1 — — arcsm 



X 



x=\J max{0,a- } 



4 



(50) 



= ^/2/i c /?ia 



where we have defined a± = 2(e ± fi c )/hu>, and m± = yo± T 

To obtain an analytic expression for the weighted density of states p w (e) we note that 



= 2^#% c (t) - tU(r)] - 

at at L ^ cV ; cxtv n dt 

We therefore find that the weighted density of states is given by 

4 



(51) 



Pw(e) 



dfi c 

dt 
dt 



TxhuO JU< 



<hn, 2 (^g - >r 



2(e - /i c 



+ y 2 



4 



/-(e)-p(el 



(52) 



The third important quantity that can be calculated analytically in the Thomas-Fermi 
approximation arises in the ergodic projection of C12. With the variable change in Eq. fl31p, 
and noting that U m i n (t) = n c {t) in the Thomas-Fermi approximation, Eq. ( j30l ) can be written 

as 
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Ji2(ei) = ^3^7 J de 2 J de 3 J de 4 [5(e x - e 2 ) - 8(e x - e 3 ) - 5(e x - e 4 )] 

X [(1 + 92)9394 - 92(1 + 03)(1 + </4)] 

x / drn c (r,t)S , (p 2 ,p3,p 4 )5(e2 - e 3 - e 4 ) 

JU<emir, 



(53) 



where Pi = y2m(ej — U). It is apparent that the integrand is symmetric in the variables e 3 
and e 4 . We can therefore assume, without loss of generality that e 2 > e 3 > e 4 which also 
implies p 2 > p 3 > Pa- In this situation, 5 , (p 2 ,p 3 ,p 4 ) in Eq.(EPI) reduces to 



S(jp%,Pz,p±) = - [1 - sgn(p2 - p 3 - , 



which is nonzero and equal to 1, if 



P2 < P3 + Pa ■ 



(54) 



(55) 



This restricts the spatial integration in Eq. (p3|) to the domain specified by this inequality. 
Inserting the definitions of pi and using the conservation of energy condition, e 2 = e 3 + e 4 , 
Eq. ( |55D is equivalent to 



F{U) = U 2 - -(e 3 + h)U + -e 3 e 4 > . 



(56) 



The roots of F(U) = are given by 

2 



e 3 + e 4 ) ± A/e| - e 3 e 4 + e\ 



(57) 



in terms of which F(U) = (U — U-)(U + U+). The requirement F(U) > is therefore 
satisfied for U < U- and £7 > C7+. The latter condition, however, is inconsistent with the 
constraint U < e 4 for the integral in Eq. (|53|). Because £7_ does satisfy U- < e 4 , the net 
effect of the factor 5'(p 2 ,p 3 ,p 4 ) in Eq. fl53D is to restrict the spatial integration domain to 
the domain defined by U < f7_, i.e., 

/ drn c (r,t)S(p 2 ,P3,P4)5{e 2 - e 3 - e 4 ) = / dr n c (r, t)5(e 2 - e 3 - e 4 ) . (58) 

JU<1^ JU<U- 



The remaining spatial integral in Eq. ( |58| ) can be carried out analytically for the Thomas- 
Fermi density profile. If U- > // c , we have simply 



f drn c (r,t) = N c (t) 

JU<U- 



(59) 



On the other hand, for < U- < fi c , we have 

/ _ drn c (v,t) = N c (t)\l 1- (l-—) 

JU<U- 2 \ fi c J 



3/2 



3 
2 



V 1 He, 



5/2' 



(60) 



Physically, Eqs. fl59D and ( 160|) are a consequence of the kinematical constraints for scattering 
into the condensate that appear in the original form of the collision integral in Eq. (|30D . 
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Finally, we indicate some implications of the assumption of adiabatic growth in the 
context of the Thomas-Fermi approximation. We take as an approximate solution to the 
dissipative nonlinear Schrodinger equation in Eq. (0) a condensate wavefunction of the form 

*(r,t) = vWr,t) e *<*\ (61) 



where n c (r, t) is the Thomas- Fermi density profile in Eq. (|46|). Inserting this wavefunction 
into Eq. ([/]), neglecting the kinetic energy term as in the Thomas- Fermi approximation, and 
separating real and imaginary parts of the resulting equation, we obtain the relations 

0(r,t) = -~^dtV c (t'), (62) 

and 

R{r ' t] --2^j7J)- (63) 

The fact that the phase is spatially independent implies that the superfluid velocity v c is 
zero as we have assumed in Sec. According to Eq. fl36|) , Eq. ( |63"D implies 

«9iV c 



Of 



-fdrfi c (t), (64) 

(J J 



where the integral is restricted to the region occupied by the condensate, i.e., n c (r, t) ^ 0. 
This is the same expression obtained by taking the time derivative of the integral of Eq. ( [46]) 
over all space. We therefore see that the wave function in Eq. (|61|) is an internally consistent 
solution of the dissipative nonlinear Schrodinger equation. 



VI. RESULTS 

In this Section we present the results of our calculations, which were performed for the 
situation corresponding to the MIT experiments [JHJ. These used 23 Na atoms confined in 
an axially symmetric trap with harmonic frequencies of 18.0 Hz and 82.3 Hz along and 
perpendicular to the symmetry axis, respectively. These values give an averaged frequency 
of u/2ir = 49.6 Hz, which implies that Tiuj/kB is equal to 2.4 nK. The s-wave scattering 
length a is 2.75 nm. 

To begin, we provide a few of the numerical details. We used a discretized energy mesh 
consisting of equally spaced points in the range < e < e max . The value of the temperature 
used in the simulations is typically of the order of 1 fiK, which requires a maximum energy 
range of about e max ~ 2500 — 3000 fuD in order to ensure that p(e)g(e) is sufficiently small 
at the end of the range. In evaluating the collision integrals I22 in Eq. (26) and in 
Eq. (30), the delta functions were used to perform some of the integrations analytically. The 
remaining integrals were then evaluated numerically using a simple trapezoidal integration 
scheme. The main advantage of this scheme in the case of the I22 collision integral is that 
the conservation of both particle number and energy is numerically exact, which in general 
is not the case for higher order integration schemes such as Simpson's rule. This conserving 
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property is especially important in simulations of the condensate growth since a loss of 
either particles or energy due to numerical inaccuracy would lead to systematic errors in the 
final equilibrium values for various physical quantities. The situation for the lyi collision 
integral is somewhat different since neither of the integrals in Eq. (35) or (36) is zero. Thus 
the numerical results will depend on the choice of the energy mesh size, and one must 
check that the results obtained for a given simulation are insensitive to variations in this 
parameter. The checks we have performed indicate that errors in the final results coming 
from this numerical source are no larger than a few percent. Errors of this magnitude will 
not influence the general conclusions that we make. As a final point, we used the Euler 
method to propagate Eq. (32) in time. The time step was chosen to be sufficiently small, 
typically 0.5 ms, to ensure the accuracy of the time evolution. 

To start the simulation we begin with an initial nonequilibrium distribution which is 
meant to represent the conditions immediately after the rapid evaporative cooling quench 
used in the experiments. Ideally, such a quench starts with an equilibrium distribution at 
some temperature T above T c and excises all particles with energy above E cut = k B T cut . 
We model this by a truncated Bose distribution at the temperature T. Although we expect 
this initial distribution to represent the experimental situation reasonably well, there will 
no doubt be differences from the actual distributions due to the finite time taken to perform 
the quench, which allows some equilibration to occur, and the possible incomplete removal 
of all particles in the energy range of the sweep. Since the rf field is resonant only at 
certain positions in the trap, atoms of a given energy must have sufficient time to reach 
these positions in order to suffer a spin flip and thus be ejected from the trap. If this is not 
the case, the distribution in energy will also have a spatial dependence. Some indication 
that such a nonergodic state in fact occurs is provided by the observation that the thermal 
cloud starts to oscillate after the quench. However, for lack of detailed information about 
the experimental initial conditions, we shall assume an idealized truncated Bose distribution 
as our initial condition. 

To complete the specification of the initial state we must also make a choice for the 
number of atoms initially in the condensate. Of course, if this number is zero, lyi as given 
by Eq. (27) is zero since we have only included stimulated transitions into the condensate. 
In the absence of spontaneous processes there is no possibility of condensate growth. Under 
the experimental conditions of interest, however, the lowest quantum state initially already 
has a rather large thermal occupation and stimulated processes will dominate. We therefore 
choose the initial condensate number to be given by the occupation of the lowest harmonic 
oscillator state at the temperature of the truncated Bose distribution. This number is 
typically of the order of a few hundred particles. As our numerical results presented below 
will show, the growth curves are rather insensitive to this starting value as long as it is small 
compared to the final equilibrium number of condensate atoms. 

In Fig. 1 we show a sequence of growth curves which illustrate the dependence on the 
parameter T cut . In this set of simulations we assume that the temperature of the equilibrium 
Bose distribution is equal to T c = 0.765 jj,K and its chemical potential jl is equal to zero. 
Before the cut, the gas contains A = 40 x 10 6 thermal atoms and the number of condensate 
atoms is given by N c = [exp(3/37kZ>/2) — = 214. In a particular simulation, the total 
number of atoms and the average energy per atom of course depends on the depth of the 
energy cut. The growth curves are characterized by an initial stage of slow growth during 
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which the truncated Bose distribution evolves into a quasi-equilibrium distribution, a well- 
defined onset time i nset where a significant increase in the rate of growth occurs, and finally 
a relaxational stage where the condensate number approaches a final equilibrium value. As 
the cut is made deeper and deeper, this final number at first increases due to the decreasing 
total energy of the initial distribution, which results in a lower final temperature. However, 
at some point the final number of condensate atoms reaches a maximum and then decreases 
with further deepening of the cut due to the reduced total number of atoms in the initial 
distribution. To distinguish this behavior the growth curves are shown as solid lines when 
the final number is increasing with decreasing T cut , and conversely, by dashed lines when 
the final number is decreasing. 

Although all the growth curves in Fig. [I] are qualitatively similar, it is clear that there 
are important differences in detail. For the curves with an increasing equilibrium number 
of condensate particles, i.e., the solid curves, both the onset time and subsequent relaxation 
time are seen to decrease with decreasing T cut . However, for the curves with an decreasing 
equilibrium number of condensate particles, i.e., the dashed curves, the dependence of both 
of these times on further decreases in T cut is much weaker, and they appear to approach 
limiting values. In order to quantify this behavior, it is convenient to fit the relaxational 
part of the theoretical growth curves to a simple exponential relaxation 



where N^ q , 7 and t nset are fitting parameters. This functional form is found to provide a 
very good fit to this part of the theoretical curves. Fig. 2 summarizes the results for the 
onset time, t nset, and exponential relaxation rate, 7, for the particular simulations presented 
in Fig. 1. The onset time decreases from about 100 ms to 20 ms as T cut /T c is reduced from 
5 to 0.5. At the same time, the relaxation rate increases from about 6 s _1 to 12 s~ 4 . 

We have also looked at the dependence of the growth curves on the other parameters 
that appear in the theory. In Fig. 3 we show the growth curves for a range of initial 
temperatures. Prior to the quench, these initial temperatures are larger than T c , and in 
each case the chemical potential is adjusted to provide again a total of 40 x 10 6 atoms in 
the thermal cloud. The energy cut and initial number of condensate atoms were taken to 
be T cut /T c = 2.5 and iV c (0) = 214, respectively, and were the same for all the runs. Not 
surprisingly, we find that the final equilibrium condensate number decreases with increasing 
initial temperature as a result of the larger average energy per atom. This of course also 
leads to a higher final equilibrium temperature. However what is somewhat unexpected is 
the very rapid increase of the onset time as the initial temperature is increased. In Fig. 4(a) 
we show that a 30% variation in T/T c gives rise to more than a ten-fold variation in t onS et, and 
that these values are typically much larger than those found using an initial temperature 
of T = T c . In addition, fig. 4(b) shows that the relaxation rate tends to decrease with 
increasing T/T c and is comparable to the values given in Fig 2. 

In Fig. 5 we show the variation of the growth curves with the initial number of condensate 
atoms. In this case, the initial nonequilibrium distribution is held fixed, corresponding to a 
Bose distribution with iV = 40 x 10 6 , T = T c and T cut /T c = 2.5. The growth curve is rather 
insensitive to the condensate number in the range 10 2 < N c < 10 4 , but then shows a much 
stronger dependence in the range 10 4 < N c < 10 6 . At the higher end of this range, the initial 
number is already visible on the graph and by N c = 10 6 there is no longer a meaningful 




(65) 
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onset time. This would correspond to a situation in which a significant condensate fraction 
has already formed by the time the quench is completed. This kind of behavior is indeed 
also seen experimentally under certain conditions. 

In order to explain some of these results it is necessary to examine the time evolution 
of the distribution function g(e,t). In Fig. 6 we show ln(g) vs. e for various times after the 
quench. At early times the distribution function is equilibrated by the scattering of thermal 
atoms into states above the E cut which are initially depleted. To conserve energy, the mean 
energy of the atoms below E cut must decrease. In fact, the population of the low energy states 
increases significantly before the onset of rapid condensate growth. This is shown in Fig. 7 
where g(e,t) is plotted as a function of time for some specific energy values. We see that 
g(e, t) at first increases rapidly, reaches a maximum at a time very close to the onset time and 
then relaxes towards its final equilibrium value of (e^ 6 — l) -1 . This behaviour is typical of 
all situations in which the growth of a condensate is observed. This strong correlation of the 
peak position in Fig. 7 with the onset time suggests that condensate formation is triggered 
by an enhanced low-energy population. Before the onset time, we find numerically that g(e) 
behaves approximately as (e)~ L63 , which is a stronger singularity than that exhibited by an 
equilibrium Bose distribution with zero chemical potential, and agrees within our numerical 
accuracy with the (e)~ 5 / 3 dependence predicted by Svistunov Regardless of the precise 



exponent, it seems that a 'super-critical' behavior of the distribution function is a precursor 



to condensate formation 24 



A useful way to characterize the time evolution of g(e, t) is to express it locally as a Bose 
distribution 

<KM) = ( n l ; , (66) 

exp(pe — //) — 1 

where the two parameters (3 and p, are defined by fitting this expression to the value of 
the distribution function and its energy derivative. Although the parameters are treated 
locally as constants in this procedure, they nevertheless depend parametrically on the energy 
variable e. The local temperature and chemical potential parameters defined in this way are 
shown in Fig. 8(a) and (b) at time intervals of 0.05 s for a situation in which the quenched 
thermal cloud equilibrates to a final temperature above T c . Both parameters are seen to 
be strongly energy dependent at early times but evolve towards energy-independent values 
by the end of the simulation. The negative equilibrium value of the chemical potential 
corresponds to an uncondensed thermal cloud at a temperature of about 1.92 /iK. 

A situation in which the quench leads to the formation of a condensate is illustrated in 
Figs. 9(a) and (b). The parameters are plotted at 0.25 s intervals during the relaxational 
stage of the growth curve beyond the onset time. At low energies, the local temperature 
lies above the final equilibrium value which reflects the higher temperature of the initial 
Bose distribution. However at higher energies, the local temperature is lower than the final 
temperature since the gas in this energy range is effectively colder as a result of the quench. 
Fig. 9(b) shows the corresponding variation of the chemical potential. As a result of the 
formation of the condensate, the chemical potential at low energies is pinned to zero and 
then increases at higher energies. The deviations of both the local temperature and chemical 
potential from their final equilibrium values are seen to relax to zero on a time scale which 
is comparable with the relaxational stage of the condensate growth. This relaxation rate 



19 



can therefore be attributed to the relatively slow equilibration of the local temperature and 
chemical potential of the thermal cloud. 

We finally turn to a comparison with experiment. This is shown in Fig. 10 for the 
particular case in which the starting number of noncondensed atoms is 40 x 10 6 , as in 
the simulations discussed above, but with the initial number of condensate atoms set to 
N c (0) = 10 4 . In the particular experimental run starting with this total number of atoms 
before the rf quench, the condensate number is found to relax to a final number of 1.2 x 10 6 
atoms. According to Fig. 1, there are two values of T cut which will lead to this final number 
of condensate atoms, T cut /T c = 0.6 and T cut /T c = 5.7. The results for the deeper cut of 
T cut /T c = 0.6 are shown as curve (b) and are seen to be in very good agreement with 
the experimental results. However, we cannot claim good agreement overall since the total 
number of atoms after the quench is only 2.5 x 10 6 as compared to the experimental number 
of about 16.0 x 10 6 atoms. For the shallower cut of T cut /T c = 5.7 shown as curve (c), the 
agreement between the theoretical and experimental growth curves is clearly worse in that 
the theoretical growth rate is too small. Moreover, the total number of atoms remaining in 
the trap is 37.5 x 10 6 which is too large by roughly a factor of 2. Alternatively, one can 
choose a cut which reproduces the final number of atoms in the trap. In our simulations, 
this requires a cut of T cut /T c = 1.9. Although the initial growth rate agrees with experiment 
in this case, the final equilibrium number of condensate atoms is 4.5 x 10 6 , which is too 
large by almost a factor of 4. This final number could be improved by elevating the starting 
temperature (recall that these simulations used T = T c = 0.765 /iK), however as Fig. 3 
shows, achieving a four-fold reduction in the equilibrium condensate number would increase 
the onset time well beyond the experimental value. It therefore appears that the present 
simulations cannot reproduce all aspects of the experiments simultaneously. 

Fig. 10 also shows a theoretical growth curve for the same initial conditions as for curve 
(b), but with mean- field interactions between the condensate and thermal cloud turned 
off. To elaborate, the potential acting on the thermal cloud is simply the time-independent 
trapping potential, and the condensate is taken to have essentially a delta-function spatial 
distribution at zero energy. In this case, the integral of n c (r, t) in Eq. (30) is replaced by 
N c (t). It can be seen that the qualitative behavior is very similar to the fully interacting 
simulation, but that the equilibrium number of condensate atoms is increased considerably, 
as expected. 

Fig. 11 provides a comparison with another set of experimental results. In this case 
the initial number of atoms before the quench is not known and was therefore taken to be 
60 x 10 6 in order to optimize agreement with experiment. Furthermore, the energy cut was 
chosen as T cut /T c = 2.5. This leads to a final number of 7.3 x 10 6 condensate atoms in the 
trap, which is approximately the same number as found in the experiment, N c = 7.2 x 10 6 . 
Although this simulation achieves good agreement between theory and experiment for the 
condensate growth curve, there are too many unknown variables, including the final number 
of atoms in the trap, to know whether or not theory is reproducing experiment. For this 
reason, the results in Fig. 11 should simply be viewed as a possible fit to the experimental 
data. 
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VII. DISCUSSION AND OUTLOOK 



Our main objective has been to obtain a realistic description of condensate growth which 
takes into account the effects of mean-field interactions. Within the ergodic approximation 
for the noncondensed atoms, and the adiabatic approximation for the condensate, the ki- 
netic equation we obtain is given by Eq. (|19"D, and we have used this equation to perform 



simulations of condensate growth. In agreement with earlier work [|18,[21|, we find that the 



growth curves have a well-defined onset time, after which an exponential relaxation towards 
equilibrium takes place. Detailed comparison with the results reported in Ref. ]10| shows 
that certain parameters can be tuned in order to achieve agreement with the experimental 
growth curves. However it seems impossible with the present simulations to reproduce the 
overall equilibrium state of the trapped gas. 

If we attribute the existing discrepancies to theory we must at some point reexamine the 
two major assumptions made in this work, namely the adiabatic growth of the condensate 
and the ergodic evolution of the thermal cloud. The adiabatic assumption neglects the 
dynamics of the condensate, specifically the possibility that collective oscillations are excited 
during the growth process. Whether or not this has any important effect on the rate at which 
atoms are exchanged between the condensate and thermal cloud in not known and should 
be investigated. In the same vein, oscillations of the thermal cloud seen in the experiments 
clearly indicate the non-ergodic state of the gas which in principle might be important in 
determining the time scale of equilibration. However, to answer this question requires a 
solution of the full quantum Boltzmann equation which seems out of reach at the moment. 
One cannot of course discount the possibility that there are uncertainties in the experimental 
results themselves. Further experimental work is needed to confirm the earlier results and to 
explore in more detail the dependences on various parameters such as the initial temperature 
of the cloud and the depth of the rf cut. 

After completion of this work, a preprint by Davis, Gardiner, and Ballagh appeared f25 



which is a continuation of a series of papers by Gardiner et al. It also addresses the issue 
of mean-field interactions as affecting the density of states, and improves on the authors' 
earlier work by giving a more realistic description of the rf quench used in the experiments. 
Thus, although there are differences in methodology, the physical basis of their work and 
the approximations they make are essentially equivalent to ours. As confirmation of this 
equivalence, their calculations of condensate growth performed for the initial conditions of 
Figs. 10 and 11 yield results which are in quantitative agreement with ours. The situation 
considered in Fig. 10 is optimal from a theoretical point of view since the experimental 
conditions are best known in this case. Yet both sets of calculations are unable to reproduce 
the experimental results in every detail. 

One of the differences between their work and ours concerns the way that the condensate 
is treated. In our formulation, the condensate is isolated explicitly as the macroscopically 
occupied quantum state, and the remaining excited states making up the thermal cloud 
are treated semiclassically. As a result of this formulation, we have two kinds of collision 
integrals, one for thermal atoms scattering amongst each other and a second for collisions of 
thermal atoms with the condensate. In the formulation of Davis, Gardiner, and Ballagh on 
the other hand, all states including the condensate are treated equivalently and thus only 
a single collision integral enters. As a result, the effective collision cross-section involving 
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the condensate does not depend on time as it does in our formulation. A second apparent 
difference has to do with the term involving the weighted density of states p w in Eq. (fl9|). 
This term arises as a consequence of the time dependence of the mean-field interaction. 
Although Davis, Gardiner, and Ballagh also deal with a time-dependent density of states, 
the second term on the right hand side of Eq. ( p~5j ) does not appear explicitly in their 
kinetic equation. However, they account for this term by dividing phase space into energy 
bins having widths which are a function of time. A final difference involves the use of 
the Bogoliubov excitation spectrum in the calculation of their density of states, instead 
of the Hartree-Fock dispersion used here. We do not expect this to affect the condensate 
growth curves significantly. However, if quasi-particle excitations are invoked, one should 
in principle also use these states to calculate the collision integrals |14[|. It is not known at 
present what effect this might have on the collision rates for the low-lying energy levels. 

Finally, we note that the ergodic treatment of the Boltzmann equation is a powerful, al- 
beit approximate, method which would allow the study of nonequilibrium processes in other 
situations as well. Some future applications might include the nonequilibrium dynamics of 
fermion-fermion and boson-fermion mixtures. Thus far, the problem of evaporative cool- 
ing in these systems has been studied using a simplified procedure whereby the distribution 
function is assumed to be given by a cut-off equilibrium distribution function PB . A cooling 



trajectory in phase space is then generated by solving for the temperature, chemical poten- 
tial and cut-off energy at each successive time step. The accuracy of this approach could 
be checked by solving for the entire distribution function following the methods used here. 
Another interesting application would be to study a nonequilibrium steady state situation 
in which atoms are continuously fed into the trapping potential while simultaneously being 
removed by a rf-cut p7|. This would be relevant to the study of steady-state atom lasers. 
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APPENDIX A: EVALUATION OF THE WEIGHTED DENSITY OF STATES. 



In this Appendix, we summarize the steps needed to include in our calculations the 
effect of the mean-field interactions arising from the noncondensed cloud itself. Referring to 
Eq. fllBD) we see that we must evaluate dU(r,t)/dt. This quantity is given by 



dU (r, t) _ 2 I dn(r,t) + dn c (r, t) " 



dt * \ dt 

The time derivative of n(r, t) can be expressed as 
dh(r,t) d f d 3 p 



Ot 



(Al) 



dt 



dtJ (2tt^) 3 
J dep(r, e,t) 
dU(r,t) 



dU(T,t)dg(e,t) dg(e,t) 
dt de dt 
dg(e,t) 



I(r,t)^V+Jdep(r,e,t) 



dt 



where we have defined 



I(r,t) = J dep(r,e,t) 



dg(e,t) 
de 



(A2) 



(A3) 



Substituting Eq. (|A1|) into Eq. (|A2|) , the latter can be rearranged to provide an expression 
for the time rate of change of n(r, t) in terms of the time rate of change of the condensate 
density n c (r,t) and the distribution function g(e). We find 



dn(r, t) 2gl(r,t) dn c (r,t) 
di ~ 1 - 2gl(r,t) di 

Inserting this result into Eq. (|A1| ), we have 
8U(r,t) 2g dn c (r,t)dN 



de 



p(r,e,t) dg(e,t) 
l-2gl(r,t) dt 



dt 



+ 



■2g 



l-2gl(r,t) dN c dt l-2gl(r,t) 



J dep(r,e,t) 



dg(e,t) 
dt 



(A4) 



(A5) 



We have here made use of the fact that n c (r, t) depends on time parametrically through 
N c (t), so that 



<9n c (r, t) dn c (r, t) dN c dn c (r, t) dN 



dt dN c dt 

Thus, the weigthed density of states becomes 

dU(r,t) 



dN r dt 



(A6) 



p w (e,t) = J d 3 rp(r,e,t)- 



dt 



J d 3 rp(r,e,t) (- 



2g dn c (r,t)\(M 
' dt 



+ J d 3 rp(r,e,t) J de 

dN r 
A (e,t)-^ + J de'B(e,e',t) 



-2gl(v,t) dN c 
f 2gp(r,e',t) dg(e',t) 
l-2gl(r,t) dt 



dt 



(A7) 
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where 



A{e,t) = - J d 3 rp(r,e,t) 



( 



1 - 2gl(r, t) dN c 



2g dn c (r,t) 



) 



(A8) 



and 



B{e,e',t) = 2g J 



3r p(r,e,t)p(r,e',t) 
l-2gl(r,t) 



(A9) 



We recover the expression for p w (e, t) given in Eq. ( p2|) by setting the kernel B equal to zero 
and neglecting I in the expression for A. It can be seen that including the mean-field of the 
noncondensate complicates the calculations considerably, but all quantities can in principle 
be calculated explicitly if these refinements are desired. However, as discussed in Sec. 0, we 
do not expect these effects to be quantitatively important. 
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FIGURES 



FIG. 1. Growth curves for different initial energy cutoffs. As discussed in Sec. |V|, the initial 
conditions are defined by fixing the temperature T = T c = 0.765 fiK and the chemical potential 
jl = of the distribution function before it is truncated. This gives N = 40 x 10 6 noncondensate 
atoms. The number of atoms initially in the condensate is chosen to be N c = 214. The solid curves 
in order of increasing saturation values correspond to T cut /T c = 5.5, 5.0, 4.5, 4.0, 3.5, 3.0 and 2.5. 
The dashed curves in order of decreasing saturation values correspond to T cut /T c = 2.0, 1.5, 1.0 
and 0.5. 



FIG. 2. The onset time (a) and relaxation rate (b) for the growth curves in Fig. [|, as 



determined by using the fitting function in Eq. (35). 



FIG. 3. Growth curves for different initial temperatures T. In order of decreasing equilibrium 
number of condensate atoms, T/T c = 1, 1.05, 1.1, 1.15, 1.2, 1.25 and 1.3. The initial conditions 
are defined by fixing the number of noncondensed particles to N = 40 x 10 6 , and the number of 
condensed atoms to N c = 214, as in Fig. @. The cutoff is now kept fixed at T cut /T c = 2.5. The 
chemical potential is less than zero, and adjusted to keep the number of noncondensed particles 
fixed. 



FIG. 4. The onset time (a) and relaxation rate (b) for the growth curves in Fig. |3j, using the 
same fitting procedure as in Fig. |2| 



FIG. 5. Growth curves for different initial number of condensed particles. The other 
parameters defining the initial conditions are the same as in Fig. |], and the cutoff is held 
fixed at T cu t/T c = 2.5. In order of increasing saturation values, the curves correspond to 
JV c (0) = 10 2 , 10 3 , 10 4 , 10 5 and 10 6 . 



FIG. 6. Plot of the logarithm of the distribution function g(e,t) for the curve with 
T cu t/T c = 2.5 from Fig. at time intervals At = 0.02 s, starting from t = 0.02 s. Each curve is 
shifted up by one unit with respect to the previous one for clarity. 



FIG. 7. Plot of the product of the density of states and the distribution function at energies of 
30, 60, 120 and 240 hCu, for the curve with T cu t/T c = 2.5 in Fig. |]. The peaks occur in the vicinity 
of the onset time i Q nset- 
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FIG. 8. Equilibration of the local temperature and chemical potential for a situation in which 
the final equilibrium temperature of the gas is above the critical temperature. Panel (a) gives the 
local temperature as a function of energy for a sequence of times during the equilibration process. 
In equilibrium, both the temperature and chemical potential are independent of energy. Panel 
(b) gives the corresponding variation of the local chemical potential. The initial conditions before 
the distribution is truncated are defined by a temperature T = 2fiK, and a chemical potential 
\x = —200. The cutoff is at T cnt /T c = 2.5. 

FIG. 9. As in Fig. [8|, but for a situation in which the final equilibrium temperature is below the 
critical temperature. The equilibration of the local temperature and chemical potential corresponds 
to the T cut /T c = 2.5 growth curve in Fig. [[]. 

FIG. 10. Theoretical growth curves for the initial conditions of Fig. [j], but with N c (0) = 10 4 
and for various energy cuts: (a) T cut /T c = 1.9, (b) T cut /T c = 0.6, (c) T cut /T c = 5.7. The exper- 
imental points are taken from Fig. || of Ref. [10]. The dashed line shows the theoretical growth 
curve for the conditions of case (b) but with mean field interactions turned off. 

FIG. 11. Theoretical growth curve for initial conditions given by N(0) = 60 x 10 6 , 
T = T c = 0.876 fj,K, T cnt /T c = 2.5 and iV c (0) = 50 x 10 4 . The experimental points are taken 
from Fig. | of Ref. [10]. 
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